(*---------------------------------------- Inputs --------------------------------------*)
(*We input an alphabet*)

Alf={1,2,3,5};
(*We input a size of matrix for the interpolation*)
m=20;
(*We imput the number of decimal places we will wotk to*)
dec=50;
(*We input an initial interval *)
t0=7/10;
t1=8/10;

(*We input a scale at which we want the programme to end*)
eps= 10^(-11);

(*----------------------------------- Setting up the interpolation machinery ----------------------------*)

(*We construct Lagrange interpolation functions ff[j,x] and the set of Chebychev points tab*)
x[k_]:= N[Cos[Pi*(2*l-1)/(2*m)],dec]/. l-> k
y[k_]:= (x[k]+1)/2 /. l-> k
tab=Table[y[k],{k,1,m}];
S[n_]:= Drop[tab,{n}]
ff[j_,y_]:= Product[(x-S[j][[k]])/(tab[[j]]-S[j][[k]]), {k,1,m-1}]/. x-> y

(*-----------------------------------Setting up an iterative loop -----------------------------------*)
n=1
While[t1 - t0 > eps, 
Print[n "th interval"]
Print[N[t0,dec]];
Print[N[t1,dec]];
(* t2 is the midpoint of t1 and t0 and f0is the interpolation test function for t2 *)
t2=(t0+t1)/2;
L[i_, y_]:=  Sum[(Alf[[j]]+x)^(-2*t2)*ff[i, 1/(Alf[[j]]+x)],{j,1,Length[Alf]}]/.x->y;
M= Table[L[i, tab[[j]]], {i,1,m},{j,1,m} ];
MT=Transpose[M];
weights= Re[Flatten[Eigenvectors[MT,1]]];
f0[y_]:= Sum[weights[[i]]*ff[i,x], {i,1,m}]/. x->y;
(* Lf0[x] is the image of f0 under the transfer operator and m0, m1 are the min and max of the quotient Lf0/ff*)
Lf0[y_]:= Sum[ (Alf[[j]]+x)^(-2*t2)*f0[ 1/(Alf[[j]]+x)],{j,1,Length[Alf]}]/.x->y;
m0=MinValue[{Lf0[x]/f0[x],0<x<1}, {x}];
m1=MaxValue[{Lf0[x]/f0[x],0<x<1}, {x}];
(*Print[m0];
Print[m1];*)
(* If m0 > 1 then we replace the lower bound t0 is replaced by t2.  Otherwise we replace the upper bound by t1 *)
If[m0 > 1,t0=t2,t1=t2];
n=n+1;
(* Now we go back to the start of the loop*)
]
(*The output is a list of shrinking intervals each containing the dimension*)