function M=qFstochastic(A,B,scope,rpts,siz,scal,opt,sav,tab,M)

%qFstochastic(2,2,10,10,1100,1,0,2); %Apollonian gasket
%qFstochastic(2,2,10,10,1100,2,1,2); %As above using the Jorgensen parametrization
%M=qFstochastic(1.991+0.05i,1.93-0.05i,120,1,1100,1,2,1,2.2-1.8i,M);
%Non-parabolic commutators using 3 params
%qFstochastic(2,2,10,10,1100,1,2,2,2.1-2i);
%Recalculating using the existing matrix 
%M=qFstochastic(1.91+0.05*1i,1.91+0.05*1i,10,1,1100,1/28,0,2,2-2i,M);

%The program is a messy attempt to generate a limit set from random
%iterations of the generators using commutators and related starting points
%Its value lies in its ability to form a very approximate representation
%of just about any limit set at all including chaoic and non-free examples
%It initially starts from 4 fixed points of the limit set
%plus 30 points selected from the period 4 orbit abinv(a)inv(b) repeating
%which helps to picture the Jordan curve features of quasi-Fuchsian sets

%A,B are two traces used to generate a conjugacy class of matrices a,b
%scope runs 10*scope^2 pseudorandom sequences scope times
%rpts adds multiple repeats of the random process
%siz (default=1100) is the image size 
%scal (default=1) scales the image
%opt=0 gives standard view, 1 Jorgensen view 
% 3 gives  Maskit parametrization. Requires trace of b to be 2
% 2 gives the double alarm 3 parameters in normal view r
%add 4 to each for the positive Grandma root
% Requires tab to be set e.g. using nonfree()
%to get the opposite grandma root add 4 to opt
%sav=1,3 saves the matrix, and a png file sav>=2 does clf first
%tab is an entry for the trace of a*b for option 3
%M is a re-input of the output matrix to add further points
%If using M as an input parameter you need to include tab=2-2i 
%in the previous parameter for a normal 2-parameter example 
%as the default trace for parabolic commutators.

if nargin<9
    tab=2-2i;
end
if nargin<5
    siz=600;
end
if nargin<6
    scal=1;
end;
if nargin<7
    opt=0;
end
if nargin<8
    sav=0;
end
if nargin<10
    M=ones(siz,siz);
end

if sav>=2
    clf;
end
tic; 
hold on
if opt>=4
    opt2=1;
else
    opt2=0;
end
opt=mod(opt,4);
[a,b]=grandma(A,B,opt2,opt,tab);
z=zeros(1,4);
j=1;
for i=1:4
    zh=fixp(a,b,i);
    if ~isnan(zh)
        z(j)=zh;
        j=j+1;
    end
end
z2=testq6(a,b);
z=[z(1:j-1) z2];
M=testq2b2(a,b,opt,10*scope,rpts,z,[-1.1 1.1 -1.1 1.1],M,siz,scal);
toc;
s=size(M);
tit=[num2str(A),', ',num2str(B),', ',num2str(scope),', ',num2str(rpts)];
tit2=[num2str(A),'-',num2str(B),'-',num2str(scope),'-',num2str(rpts)];
z=zeros(1,s(1)*s(2));
k=1;
sc11=1;
for i=1:s(1)
    for j=1:s(2)
        if M(i,j)==0
            z(k)=((i-s(1)/2)*2/s(1)+1i*(j-s(2)/2)*2/s(2))*sc11;
            k=k+1;
        end
    end
end
plot(z(1:k-1),'.','MarkerSize',1,'Color','k');
axis([-sc11 sc11 -sc11 sc11]);
text(-1,1.05,tit);
if mod(sav,2)==1
%imwrite(255*M',[tit,'.tiff']);
imwrite(255*M',[tit,'.png']);
save([tit2,'.mat'],'M');
%saveas(gcf,[tit,'.eps'],'eps');
end
        
function M=testq2b2(a,b,opt,its,repts,z,D,M,siz,scal)
if nargin<7
    D=[-2 2 0 2];
end
if nargin<4
    repts=1000;
end
if nargin<3
    opt=0;
end
for i=1:repts
    %tic;
    rng('shuffle') 
    M=quasi2f(a,b,its^2,its,z,D,opt,M,siz,scal);
    %toc;
    %fprintf('%d\n',i);
end

function M=quasi2f(a,b,its,repts,z,A,opt,M,siz,scal)
siz=siz/2;
sizs=10/11*siz;
if nargin<3
its=1000000;
end
if nargin<4
repts=10;
end
if nargin<7
    opt=0;
end
rh=zeros(1,4);
L=length(z);
Z=z;
for i=1:L
    zs=scal*z;
    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))))=0;
           if ~opt
              M(siz+floor(sizs*real(-zs(i))),siz+floor(sizs*imag(-zs(i))))=0;
           end
        end
    end
end
stps=its/repts;
for k=0:repts-1
    for j=1+k*stps:stps+k*stps
        if j==1
            z=Z;
        else
            z=zh;
        end
        rr=floor(4*rand());
        kh=4;
        if j>1
            while rr==mod(rh(4)+2,4)
                rr=floor((4+kh)*rand());
            end
            if rr>3
                rr=rh(4);
            end
        end
        switch rr      
             case 0
                z=(a(1,1).*z+a(1,2))./(a(2,1).*z+a(2,2));
             case 2
                z=(-a(2,2).*z+a(1,2))./(a(2,1).*z-a(1,1));
             case 1
                z=(b(1,1).*z+b(1,2))./(b(2,1).*z+b(2,2));
             case 3
                z=(-b(2,2).*z+b(1,2))./(b(2,1).*z-b(1,1));
        end
        zh=z;
        rh(1:3)=rh(2:4);
        rh(4)=rr;
        for i=1:L
            zs=scal*z;
            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))))=0;
                   if ~opt
                      M(siz+floor(sizs*real(-zs(i))),siz+floor(sizs*imag(-zs(i))))=0;
                   end
                end
            end
        end
    end
end

function z=testq6(a,b)
[x,y]=testq4(a,b,'0123','w',1000,20,0);
z=nexus(x,y,2);
[x,y]=testq4(a,b,'3210','w',1000,20,0);
z2=nexus(x,y,0);
z=[z z2];

function zzzz=nexus(x,y,m)
z=x+1i*y;
j=1;
for i=1:length(z)
    if mod(i,4)==m
        zz(j)=z(i);
        j=j+1;
    end
end
zzz=logspace(0,log(length(zz))/log(10),15);
for i=1:length(zzz)
    zzzz=zz(floor(zzz));
end

function [x,y]=testq4(a,b,s,c,its,siz,pl)
if nargin<7
    pl=0;
end
x=zeros(1,its);
y=x;
z=fixp(a,b,0);
x(1)=real(z);
y(1)=imag(z);
for j=1:its
    z=x(j)+1i*y(j);
    rp=s(mod(j,length(s))+1);
    if strcmp(rp,'0')
        rr=0;
    end
    if strcmp(rp,'1')
        rr=1;
    end
    if strcmp(rp,'2')
        rr=2;
    end
    if strcmp(rp,'3')
        rr=3;
    end
    switch rr      
         case 0
            z=(a(1,1)*z+a(1,2))/(a(2,1)*z+a(2,2));
         case 2
            z=(-a(2,2)*z+a(1,2))/(a(2,1)*z-a(1,1));
         case 1
            z=(b(1,1)*z+b(1,2))/(b(2,1)*z+b(2,2));
         case 3
            z=(-b(2,2)*z+b(1,2))/(b(2,1)*z-b(1,1));
    end
    x(j+1)=real(z);
    y(j+1)=imag(z);

end
if pl
plot(z,'.','MarkerSize',siz,'Color',c);
plot(z(1),'.','MarkerSize',siz,'Color','g');
plot(z(end),'.','MarkerSize',siz,'Color','r');

end

function z=fixp(a,b,ri)
        if ri==0
            ri=floor(4*rand()); %find a fixed point from one of the transforms
        end
        %ri=5;
        switch ri
            case {0,4}
                z=(-(a(2,2)-a(1,1))+sqrt((a(2,2)-a(1,1))^2+4*a(1,2)*a(2,1)))/(2*a(2,1));
            case 1
                z=(-(a(2,2)-a(1,1))-sqrt((a(2,2)-a(1,1))^2+4*a(1,2)*a(2,1)))/(2*a(2,1));
            case 2
                z=(-(b(2,2)-b(1,1))+sqrt((b(2,2)-b(1,1))^2+4*b(1,2)*b(2,1)))/(2*b(2,1));
            case 3
                z=(-(b(2,2)-b(1,1))-sqrt((b(2,2)-b(1,1))^2+4*b(1,2)*b(2,1)))/(2*b(2,1));
            case 5
                c=a*b*a^-1*b^-1;
                z=(-(c(2,2)-c(1,1))-sqrt((c(2,2)-c(1,1))^2+4*c(1,2)*c(2,1)))/(2*c(2,1));
        end

function [a,b]=grandma(ta,tb,posroot,opt,tab)
if nargin<5
    tab=2-2i;
end
if nargin<4
    opt=0;
end
if nargin<3
    posroot=false;
end
if opt~=2
    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 || opt==3
    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
    a=[ta-tb/tab ta/tab/tab;ta tb/tab];
    b=[tb-ta/tab -tb/tab/tab;-tb ta/tab];
    end
    if opt==3
        c=-1i/(2^(1/2));
        p=(1-c^2)/(1i*c);
        q=-c;
        r=c;
        s=1i*c;
        C=[p q;r s];
        a=C*a*C^-1;
        b=C*b*C^-1;
    end
 else
    tC=ta^2+tb^2+tab^2-ta*tb*tab-2;
    Q=sqrt(2-tC);
    if abs(tC+1i*Q*sqrt(tC+2))>=2
        R=sqrt(tC+2);
    else
        R=-sqrt(tC+2);
    end
    zo=(tab-2)*(tb-R)/(tb*tab-2*ta+1i*Q*tab);
    a=[ta/2 (ta*tab-2*tb+2i*Q)/((2*tab+4)*zo);(ta*tab-2*tb-2i*Q)*zo/(2*tab-4) ta/2];
    b=[(tb-1i*Q)/2 (tb*tab-2*ta-1i*Q*tab)/((2*tab+4)*zo);(tb*tab-2*ta+1i*Q*tab)*zo/(2*tab-4) (tb+1i*Q)/2];
end