(* 
	(c) 2010 Alex Kontorovich

Algorithm due to Jenkinson-Pollicott

 *)

(* Parameters *)

alph={1,2}; (* this is the alphabet of allowed letters in the continued fraction expansion *)


(* CODE *)
w[i_]:=Module[
{j,itemp,n,ans},
itemp=i;
ans=1;
For[j=0,j<Length[i],j++,
n=itemp[[Length[i]]];
itemp=Drop[itemp,-1];
PrependTo[itemp,n];
ans=ans FromContinuedFraction[{0,itemp}];
];
ans=Simplify[ans];
ans
];

CC[n_,s_]:=Module[{tups,ws,ans},
tups=Tuples[alph,n];
ws=w/@tups;
ans=2/n Total[ws^(2s)/(1-ws^2)];
ans
];


(* make b's *)
b2=-CC[2,s];
b4=-CC[4,s]+1/2 CC[2,s]^2;
b6=-CC[6,s] + CC[2,s] CC[4,s] -1/6 CC[2,s]^3;

(* Test their accuracy *)
Plot[1+b2+b4+b6,{s,0,1}]
FindRoot[1 +b2,{s,1/2}]
FindRoot[1+b2+b4,{s,1/2}]
FindRoot[1+b2+b4+b6,{s,1/2}]

(* Looks good *)

(* Now test alphabet {2,4,6} *)

alph={2,4,6,8,10};


(* make b's *)
b2=-CC[2,s];
b4=-CC[4,s]+1/2 CC[2,s]^2;
b6=-CC[6,s] + CC[2,s] CC[4,s] -1/6 CC[2,s]^3;

(* Test their accuracy; this takes a minute *)
Plot[1+b2+b4,{s,0,1}]
FindRoot[1 +b2,{s,1/2}]
FindRoot[1+b2+b4,{s,1/2}]
FindRoot[1+b2+b4+b6,{s,1/2}]
Beep[]